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It is known that t 



dispersing obstacles 

a 



re non-equilibrium version of the Lorentz gas (a billiard with 
l|, electric field and Gaussian thermostat) is hyperbolic if the 



field is small [2] . Differently the hyperbolicity of the non-equilibrium Ehrenfest gas 
constitutes an open problem, since its obstacles are rhombi and the techniques so far 
developed rely on the dispersing nature of the obstacles . We have developed 

analytical and numerical investigations which support the idea that this model of 
transport of matter has both chaotic (positive Lyapunov exponent) and non-chaotic 
steady states with a quite peculiar sensitive dependence on the field and on the 
geometry, not observed before. The associated transport behaviour is correspond- 
ingly highly irregular, with features whose understanding is of both theoretical and 
technological interest. 
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I. INTRODUCTION 



The Ehrenfest model of diffusion (named after the Austrian Dutch physicists Paul and 
Tatiana Ehrenfest) was proposed in the early 1900s in order to illuminate the statistical 



interpretation of the second law of thermodynamics anc 



Boltzmann equation. In the Ehrenfest wind-tree model [4j, the point-like ("wind") particles 



to study the applicability of the 



aced fixed square scatterers ("tree"). 
5| to prove that microscopic chaos is not nec- 

has been considered in 



move on a plane and collide with randomly p 
This model has been recently reconsidered in 

essary for Brownian motion. A one-dimensional version of this mode 
|6j to investigate the origin of diffusion in non chaotic systems. In 
two sufficient ingredients for diffusive behavior in one-dimensional, non-chaotic systems: a) a 
finite-size, algebraic instability mechanism, and b) a mechanism that suppresses periodic or- 



6|] the authors identify 



bits. A nonequilibrium modification of the model, with regularly placed scatterers, has been 
proposed in [71] to test the applicability of the so called fluctuation relation (5,1^, 10, 111 ]) 
to non chaotic systems. This modified model was chosen under the assumption that, at small 
and vanishing fields at least, it must be non-chaotic, since collisions with flat boundaries do 
not lead to exponential separation of nearby phase space trajectories. 

However the question of whether such a model can have positive Lyapunov exponents, as 
functions of the field, is open. Indeed, the techniques so far developed, e.g. by Chernov and 
Wojtkowski , rely on the dispersing nature of the billiard obstacles. 
In this paper, the dynamical properties of the nonequilibrium version of the Ehrenfest gas 
are considered as functions of the field and of the parameters which determine the billiard 
table. Numerical tests are performed to find chaotic attractors and to compute the Lya- 
punov exponents. The construction of a sort of bifurcation diagram of the attractor as a 
function of the electric field and of the geometry is attempted. The result turns out to be 
quite peculiar: chaotic regimes with an extremely sensitive dependence on the parameters 
appear possible, although not easy to establish rigorously. If this model can be taken to 
approximate the transport of matter in microporous membranes, our results confirm the 
sensitive dependence of microporous transport on all relevant parameters observed e.g. in 
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131 ]. Indeed, the current of the nonequilibrium Eherenfest gas is proportional to the 
sum of the Lyapunov exponents (cf. Eq. fllll) below), which varies with the parameters as 
irregularly as the attractors do. 
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II. THE MODEL 

The billiard table consists of rhombi of side length I with distances along the x and 
y directions between the centers of two nearest neighbouring rhombi given by xl and yi, 
respectively. The centers of the rhombi are fixed on a triangular lattice in a plane and have 
coordinates 




= m c \i + n c \ 2 , m c , n c G Z 



where li = (xl,0) and I2 = (0, yz) are the lattice vectors. If all the pairs (m c , n c ) are 
selected, the billiard is invariant under the group of spatial translations generated by li and 
I2. Accordingly, the whole lattice can be mapped onto a so-called Wigner-Seitz cell, with 
periodic boundary conditions (Figured]). The elementary Wigner-Seitz cell of the triangular 
lattice is a hexagon of length side L and area 

A W s = |li x 1 2 |. 

The centers of all other cells are identified by the pairs (m c , n c ) G 7L E x 7* E or (m c , n c ) G 
Z° x Z° where 7L E = {n G Z : \n \ is even}, Z° = {n G Z : \n \ is odd}. Because of 
the bi-jective correspondence between rhombi and pairs (jn c ,n c ), one may label a generic 
rhombus by TZ mc , nc and the corresponding hexagon by H mcjnc . Further, a label can be put 
on the sides of the rhombi and of the hexagons, introducing an alphabet A = {ri, r 2 , r%, r^}, 
starting from the right vertical side and oriented clockwise, for the sides of the rhombi and 
an alphabet B = {hi, h 2 , h 3 , hi, h 5 , h 6 } for the hexagon sides, starting from the right vertical 
side and oriented clockwise. In this alphabet, the sides of a generic rhombus of the lattice 
can be labelled by a triple (TZ mcjnc ,s) with s G A, while the hexagon side with a triple 
(7~(-m c ,n c , s) = (m c ,n c ,s) with s G B. The rhombi lying in the y axes have TZ 0ii = (0,i) with 
i G Z° and the ones lying in the x axes have TZjfi = (j, 0) with j G Z E . 
The geometry of the model is determined by the side L of the hexagonal Wigner-Seitz cell, 
so that xl = V%/2 and y^ = 3L/2. Let TZq = (0, 0) be the rhombus with the center in the 
origin of the Cartesian coordinates, I its sides length, s x and s y the half length of the major 
and minor diagonals respectively, so that I = */ ' s 2 x + s^. To prevent the overlap of rhombi, 
the side length of the rhombus inside one hexagon has to verify 
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~ ~ 2 
V7 

which implies < s x < xj s . The case with I = — L corresponds to a billiard table which 
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id 



was recently considered in 
Take / E ^0, -— Lj, hence s x < x/. The horizon of the billiard depends on the quantity s y 
and in particular on the difference y^ — 2s y . If s y > ?/l/2, the horizon is finite; if s y < ?/l/2, it 
is infinite. The infinite horizon case allows collision-free trajectories, parallel to the x— axis. 
When the dynamics is followed within to the Wigner-Seitz cell, the position of the point 
particle of mass M must be supplemented by the couple (m c ,n c ), in order to determine 
its actual position in the infinite plane. The space between the rhombi forms the two- 
dimensional domain, in which the particle moves with velocity v during the free flights, 
while collisions with the sides of the rhombi obey the law of elastic reflection. 
In order to drive the model out of equilibrium, an external electric field parallel to the x-axis, 
E = ex is applied. If there were no interaction with a thermal reservoir, any moving particle 
would be accelerated by the external field, on average, leading to an indefinite increase 
of energy in the system, and there would be no stationary state. Therefore, in Q the 
particle has been coupled to a Gaussian thermostat. The resulting model, with periodically 
distributed scatterers, has been called the non- equilibrium Ehrenfest gas. Its phase space 
has four coordinates (x,y,p x ,p y ) and its equations of motion are given by 

x = p x /M , p x = e - a(p)p x 

with a(p) = -ep x (1) 

y=p y /M, p y = -a{p)p y 
where e is the electric field and a the Gaussian thermostat. The quantity — div(p, q) = 
ot = ep x is known as the phase space contraction rate. Because of the Gaussian thermostat, 
Px + Py 1S a constant of motion, hence there are only three independent variables, and one 
may replace p x and p y by the angle 9 G (— tt, 7r] that p = p(cos 8, sin 9) forms with the x-axis. 
For sake of simplicity, we set M = 1 and p = 1 . 

Then, if (x t ,yt) denotes the position at time t and 6 t the velocity angle, measured with 
respect to the x-axis, the trajectory between two collisions reads jl4j ] 



%t = %o ~ 7 m 



1 i singj: 



sin ft 







yt = yo - \{e t - e t0 ) (2) 

tan y = exp(— e(t — t )) tan 
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where to is the time of the previous collision, while the collision map C is given by 



y' t = vt (3) 

where 9 t is the incidence angle, (x t , yt) is the collision point and 9 is the angle that the side of 
the rhombus makes with the x— axis. The sign ± depends on the side on which the bounce 
occurs. Hence C is piecewise linear in 9 t . 

Considering the dynamics as a geodesic flow on a Riemann manifold, the appropriate metric 



for this system is Jl5 



16, 



id 



ds 2 = e~ 2tx (dx 2 + dy 2 ) 

which implies that the quantities n y = e~ ex p y and (fi t = 9 t + ey t = 9 to + e9 are conserved. 
Also, the path length £(P ,P t ) between P = (xo,yo,9t ) and P t = (x t ,yt,9 t ) turns out to 
be 



-e ex ° sin 9 tQ \ cot 9 to - cot 9 t 



= -e 
e 



(4) 



III. PERIODIC ORBITS AND STABILITY MATRICES 

Using the symbols introduced above, a trajectory segment Qn which con- 
sists of collisions can be labelled by a finite symbolic sequence such as 
(TZ h>jl , Sx)(TZ i2jj2 , s 2 ) . . .{TZ iNtjN ,s N ), with (i c ,j c ) e Z E x Z E or (i c ,j c ) E Z° x Z° and 
s c G A for c = 1, . . . , N. Periodic trajectories will be labelled by sequences which are in- 
finitely many copies of a fundamental finite sequence. 

There are two types of periodic orbits; those that are periodic in the plane, i.e. that return 
to the initial point in the plane (they are closed: Axi = 0, Ayi = 0) and those whose peri- 
odicity relies upon the periodicity of the triangular lattice, and return to the same relative 
position in a different cell (they are open: Axi ^ or Ay^ ^ 0). 

The velocity vectors of the closed orbits with two collisions have to be orthogonal to both 
sides of the rhombi where collisions occur. This implies 9q = | + 9 and 9' Q = | — 9, where 9q 
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is the out-going velocity angle and 9' the in-coming angle, because the absolute value of the 
velocity angle decreases and preserves the sign during the free flight jl4]. Then, the closed 
period-two orbits fly between rhombuses in the same line parallel to the y-axes, and have 
period r and length £ given by |18 ] 



e 



tan ( f + | 



tanff 



C = -e~ eX0 sin6 . 
e 



Now, fix the geometry of the model through the parameters L > 0, s x > 0, s y > and take 
the initial conditions 

jx , y = ±yi, 9 

It is easy to show that the periodic orbits (R-ofi, r^)(7lo^, r 3 ) (Figure [21 left panel) and 
(TZofi, ?"3)(7^o,-2, r 4 ), with electric field 

—9 + tan 6 In cos 9 



2 



exist if and only if [18] 



29 29 
3L <e< 3L - 2s 



(5) 



(6) 



For open orbits with two collisions in the finite horizon case, simple algebra 
shows that the following symbolic representations (JZij, r 4 )(7^_i J+ i, r 2 )(7^ + 2,j, r 4 ) and 
(R-ij, Ti)(^i+3j+i, r 3 )(lZi + 2j, ri) cannot be realized. However, orbits with sym- 
bolic representation (R-ij, r 4 )(7?. i+ i ) j + i, r 3 )(7£ i+ 2j, r 4 ), and its symmetric counterpart 
(R-ij, r 3 )(7?.j + i i j_i, r 4 )(7?.j + 2j, r 3 ) (Figure [2], right panel) do exist. Their period is given by 

r = - In . 

e tan(f - 9) 

The open periodic orbits with four collisions (Figure [51 left panel) and sym- 
bolic sequence (7^o,o, r^)(7l- r 2 )(7?.i i i, r 3 )(7^o,o, r i)(^2,o, ^4), and its simmetric image 
(Ko, , r 3 )(TZ-i-i,ri)(Tli-i,r4) (fco,o, r 2 ) (7£ 2 ,o,r 3 ) do exist if [18] 



max 



f 1 . sinflni . r 1 , sin#n 1 

<^ - Sa,, -X; + - In -— - f < x < mm <^ - x t + s x + - In -— -, > 
l e sinf/nJ I. e sin n J 



(7) 



To compute the Lyapunov exponents for these orbits and any other trajectory, consider 
the stability matrix J5 for a trajectory as the product of free flight stability matrices Jf and 
collision stability matrices Jc 

n 

Js = X\Jc{i)J F (i) 



i=i 
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The number of degrees freedom for the billiard map is two and the variables that we will 
use are (x ,9 ). Thus 

/ ' ' -10 



J, 



c 



1 



96*o dxo 
dx' dx' 

The free flight matrix Jp depends on the side which the trajectory leaves and the one which 
it reaches. There are two different types of side, the ones with positive slope, of equation 
y = tan 9 x + c and the ones with negative slope, of equation y = — tan 6 x + d, where c and 
d are real numbers. 

Let us compute the free flight matrix J +i _ of a trajectory which goes from a side with 
positive slope to a side with negative slope. Let (xo,yo,9o) and (x' ,y' ,9' ) be the initial 
condition on a side with positive slope and the final condition on a side with negative slope 
respectively. By using the equations of the trajectory and by the implicit function theorem 
1 18l 1 2 31 ] we obtain the Jacobian matrix of the free flight: 

/ (tan 9 + tan 9) tan 9' 2e tan 9 tan 9' \ 



(tan 9' + tan 9) tan 9 tan 9' + tan 9 



1 tan 9n — tan 9' n tan 9' n — tan 9 



(8) 



y e tan 9 tan 9' + tan 9 tan 9' + tan 9 ) 
Similarly, the flights from a side with negative slope to a side with positive slope yield 

/ (tan 9q — tan 9) tan 9' 2e tan 9 tan 9' \ 



J_ 



(— tan^Q + tan^) tan^o — tan + tan ^ 



V 



1 tan 0q — tan # 
e tan ^ tan 9 — tan #q 



(9) 



and those from one side to a parallel one yield 

/ (tan^o T tan^) tan#Q 
(tan *5g =f tan 9) tan ^ 



J. 



±,± 



tan 9' + tan ^ 
tan 9' — tan J 

\ 







tan^n — tan^n 



(10) 



y e tan 9 tan 0q =)= tan 9 J 



If y^i and \i2 are the eigenvalues of the stability matrix J$ for a periodic orbit of period r, 

1 



the two Lyapunov exponents are Aj 



log 



r 

Ax 
r 



, z = 1,2, and one obtains 
(Ai + A 2 ) 



(11) 
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where j is the current and Ax is the corresponding displacement in the direction of the field 
jl^ . Both the Lyapunov exponents of the closed periodic orbits, with period two, vanish. 
Indeed, consider that this periodic orbit has Ax = 0, hence Ai + A2 = 0. Furthermore, the 
stability matrix of these periodic orbits, which is Js = JcJ-+JcJ+-, is given by 



/4tan 2 # 



Js 



tan 2 



(l + tan 2 0) 2 
4etan6»(tan 2 6> - 



V 



'1 + tan 2 



4etan6>(tan 2 6> - 1) \ 
(1 +tan 2 fl) 2 

4tan 2 #- (1 -tan 2 ^ 
(l + tan 2 #) 2 J 



whose determinant is 1, while its trace vanishes. This, implies that both Lyapunov exponents 
vanish. 



IV. NUMERICAL ESTIMATES OF LYAPUNOV EXPONENTS 

In this paper, a system is called chaotic when it has at least one positive Lyapunov 
exponent. We note that the boundary of our system is not defocussing and the external 
field has a focussing effect, so the overall dynamics should not be chaotic in general, although 
it is not obviuos that this is the case for all values of the electric field e. In this section we 
examine the stationary state and the Lyapunov exponents, obtained by using the algorithm 



developed by Benettin, Galgani, Giorgilli and Strelcyn |l9|], for different values of e, ranging 
from small to large fields. 

A. Chaos for large electric fields 

Numerical simulations of the model starting with random initial conditions and electric 
field in the range [0.02, 1] have been initially performed for a trajectory of length n = 10 7 
collisions. The parameters chosen for the geometry are L = 1.291, s x = 0.7573 and s y = 1.1 



which correspond to a case in which t 



re angles of the rhombi are irrational w.r.t. 71", then, 



according to a conjecture by Gutkin [20], the equilibrium version of this model, (i.e. the 
e = case), should be ergodic. For simulations of 10 7 collisions, Figure [3] shows that the 
fields which appear to lead to one positive Lyapunov exponent cover a range larger than 
that which appears to correspond to two negative exponents. However, do 10 7 collisions 
suffice for a generic trajectory to characterize the stationary state? For the cases with 
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two negative exponents the answer is affirmative, since the trajectory is clearly captured 
by an attracting periodic orbit. But the cases with one apparently positive exponent are 
not equally clear. As already noted, the doubt is that, starting from a generic initial 
condition, convergence to the steady state might be too slow to be discovered, for reasons 
which had not been investigated. Indeed, even in cases in which convergence is observed, 
the particle often appears to peregrinate in a sort of chaotic quasi steady state for very- 
long times, before eventually settling on a periodic or quasi-periodic steady state. Ref. |7J] 
suggested this might always be the case. 

Therefore, we extended the simulations of the cases with one apparently positive Lya- 
punov exponent, up to 1.5 x 10 8 collisions. The behaviour of the system still appears to be 
non trivial, in cases such as that of e = 0.374. Furthermore, plotting the last 10 4 iterates 
of a trajectory of length 5 x 10 7 collisions, we cover a large fraction of the phase space, 
which appears quite close to that covered by the last 10 4 iterates of a trajectory of 1.5 x 10 8 
collisions (Figure H]). 

The conclusion that a chaotic stationary state has been reached seems reasonable in this 
case, as the computed positive Lyapunov exponent also indicates, having apparently con- 
verged to 0.144 with three digits of accuracy, after only 10 6 collisions. 
To strengthen this results, we have looked for an unstable periodic orbit embedded in the 
attractor, and we have found one periodic orbit of period four, which apparently lies in 
the attractor and has a positive Lyapunov exponent (Figure [5j right panel). The other 
possibility is that the orbit is isolated, and is separated from the attractor by such a small 
neighborhood that is numerically impossible to resolve. 

Another interesting example is given by e = 0.5: the last 7000 points of a trajectory of 
length 5 x 10 7 collisions are compared with the last 10 4 points of the trajectory of length 
2 x 10 8 (Figure [6]). Also in this case the stationary state seems to have been reached, and 
a periodic orbit of period four with one positive Lyapunov exponent seems to be embedded 
in the attractor, similarly to the case of e = 0.374. The values of e, of the initial conditions 
and the Lyapunov exponents are reported in table I and II. 



10 



B. Chaos for small electric fields 

Numerical simulations of the model were performed starting with random initial condi- 
tions and considering the electric field in the interval [0.002,0.1], for trajectories of length 
n = 10 7 . Looking at Figure [7] (left panel) we find more cases with one apparently posi- 
tive Lyapunov exponent than with two negative exponents. Plotting the last 10 points, in 
trajectories of length 1.5 x 10 8 , for the electric fields that produced apparently positive Lya- 
punov exponents, the situation is pratically the same as for large fields: some cases with one 
exponent that appeared to be positive after 10 7 collisions eventually produce (after 1.5 x 10 8 
collisions) two negative exponents and a periodic or quasi-periodic steady state |7j] |24j ] . We 
have further increased the number of collisions up to 3 x 10 8 , but the exponent Ai remained 
positive in most of the cases (Figure L7J right panel) and an apparently chaotic attractor was 
reached. 



C. A small basin of attraction 

The behaviour illustrated above is rather peculiar and calls for some explanation. How 
can it be that a steady state is so hard to reach in so many cases? Usually, convergence to 
an attracting periodic orbit occurs rather quickly, while doubts remain in some of the cases 
we considered even after 10 8 iterations of the bounce map. Therefore, we have investigated 
in greater detail the specific example with e = 0.087, L = 1.291, s x = 0.7573 and s y = 1.1. 
One finds that the largest time dependent Lyapunov exponent, Ax, rapidly settles on a 
positive value, as if the trajectory had reached a chaotic attractor. However, for randomly 
chosen initial conditions, a striking and precise monotonic 1/n behavior sets in for Ai, after a 
critical, typically large, time N c (cf. left panel of Figure E]), as if the trajectory had eventually 
collapsed on an attracting periodic orbit. Indeed, the asymptotic value of Ai is —0.004838 



with an estimated error not larger than 10 6 [25]. The accuracy of this result is due to the 



fact that a sufficiently long simulation approaches an attracting periodic orbit to practically 
full numerical precision; the Lyapunov exponents can then be computed with analogous 
accuracy. This is particularly true in the present example, which turns out to have an 
attracting orbit made of only 19 points, whose initial condition, up to 12 digits accuracy, 
is given by x po = 0.418447478686, y po = 0.492193019206, 6 po = 0.718794450586, (cf. right 
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panel of Figure [8]). For such a small number of points, numerical errors cannot appreciably 
affect the result. In particular, no doubts remain, in this case, about the negative sign of 
both exponents, hence about the attracting nature of the orbit. 

The question now arises as to the shape and size of the basin of attraction B of the 
asymptotic periodic orbit, because a particle with random initial conditions wonders around 
almost all phase space before falling inside this set. 

Our analysis of the evolution of trajectories, with initial condition close to the attracting 
periodic orbit, shows that B is quite limited in size, and particularly hard to reach because 
it contains a very small region around the right vertex of the rhombus. Furthermore, by 
varying the initial conditions around (x po , y po , 9 po ), and computing the Lyapunov exponents, 
the irregular shape of B is evidenced by the times N c , which vary most irregularly from O(10) 
to O(10 8 ), O(10 8 ) being typical for random initial conditions. Nevertheless, the radius of B, 
i.e. the supremum distance between any two points of B, is not smaller than 10 -4 , which 
is quite small but well above the distances which can be accurtely measured with double 
precision numerical simulations. 



We conclude that B lies at the border of a chaotic repeller [2l|] , and that it is quite small 
and of irregular shape. Hence, only after a sufficiently long peregrination in phase space, 
when the transiently chaotic trajectory gets sufficiently close to B, may numerical errors let 
the trajectory jump inside B. At this stage, the sudden (exponential) convergence of the 
trajectory and the consequent 1/n behaviour of the Lyapunov exponents begin. 

The same qualitative behaviour has been observed for e =0.002, 0.003, 0.004, 0.005, 0.006, 
0.007, 0.010, 0.022, 0.040, 0.041, 0.042, 0.043, 0.050, 0.064, 0.065, 0.073, 0.087. 

V. BEHAVIOUR OF THE ATTRACTOR 

It is interesting to understand the behaviuor of the steady state with the field, i.e. to 
build a kind of bifurcation diagram, e.g. to compare with the one for the Lorentz gas, whose 
obstacles are defocussing 1141 1. Our analysis reveals substantial differences from the case of 

n 

[14], as well as from standard low dimensional dynamics. 
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A. Multi-furcation as function of the electric field 

In order to visualize the behaviour of the attractor as a function of the electric field, we 
consider a projection of the billiard map phase space: the projection onto the 9 axis, which 
shows a sort of "multi-furcation" scenario. When the field is varied, a series of dramatic 
changes in the dynamics occurs. 

For the geometry determined by L = 1.291, s x = 0.7573 and s y = 1.1, numerical simulations 
were performed for a random initial condition, ignoring the initial transient behaviour, for 
the electric fields in the range [0.01, 1] with a step Ae = 0.01. At the end of a trajectory 
of 10 7 collisions, the last 10 3 points were plotted (Figure El presents only the upper half 
projection onto the 9 axis; the other half of this diagram is trivially obtained from this, by 
reflection along the line 9 = 0, as a result of the symmetry imposed on the system by the 
external field). In this way, many electric fields are found whose dynamics sample most of 
the 9 space, while only a few fields show a periodic (or quasi periodic) steady state. The 
apparently chaotic regions and the regular regions are finely interspersed with each other, 
in a way which we have not found elsewhere in the literature. 

We have performed simulations also in the range [1, 1.3] and have only plotted the last 
10 3 points out of 10 7 collisions, as in the previous case. In this range, all attractors have 
been found to be periodic orbits, which are reached well before 10 7 collisions. In other 
words, the stationary state is rapidly reached with large electric fields as expected for the 
correspondingly high dissipations. 

How conclusive are these results? Again, when a periodic steady state is reached, the 
situation is clear, while doubts remain when the steady state looks chaotic. Therefore we 
have considered the range [0.77,0.86]: a rather wide range, apparently chaotic, but also 
characterized by high dissipation, which favours ordered dynamics. Running trajectories of 
10 8 collisions and plotting the last 10 3 points, the apparently weakly chaotic states survived, 
despite the high dissipation. This study combined with the analysis of the largest Lyapunov 
exponent makes chaos quite plausible for these fields, although we cannot exclude that the 
trajectories collapse on periodic orbits after much longer times. On the other hand, there 
is no purpose in pushing further this analysis, since it is bound to remain uncertain. Here, 
it suffices to have uncovered a rather peculiar behaviour, unexpected for simple dynamical 
systems. 
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Indeed, in our model, the passage from low-period attractors to apparently chaotic steady 
states, is rather abrupt and seems to be discontinuous. Moreover, the periodic attractors 
always seem to coexist with transiently chaotic states. Differently, standard bifurcation 
scenarios are characterized by a gradual growth of the period of the attracting orbits, the 
orbits are globally attracting and lead to fast convergence of the trajectories towards them. 
To further investigate this behaviour, we have honed the range [0.77, 0.78] with step 0.001 
simulating trajectories of 10 s collisions, finding a (quasi) periodic orbit at e = 0.771 and 
another one at e = 0.777, with apparently chaotic trajectories in the middle. Thus, we have 
further honed the range [0.770; 0.771] with a step of 0.0001, to find that a jump from (quasi) 
periodic orbits to apparent chaos happens for a field variation of just 10~ 5 . 
We have also performed simulations in the range [0.771,0.772], namely from the (quasi) 
periodic orbit at e = 0.771 to an apparently chaotic case at e = 0.772. Plotting the last 
10 3 points of trajectories of 10 8 collisions, apparent chaos and (quasi) periodic orbits appear 
again to be finely intertwined. The range [0.18,0.38] was similarly studied, and some of 
the cases that appeared to be chaotic after 10 7 collisions, turned into (quasi) periodic after 
10 8 , but not all. Therefore, the duration of the transients results exceedingly long in all 
cases, which is a manifestation of either very small basins of attraction, or of the smallness 
of the stable islands, whose size would then vary wildly with the field. Although this is 
not conclusive evidence, our analysis supports the idea that there could be a discontinuous 
transition between chaos and regular motion, in the non-equilibrium Ehrenfest gas. 

B. Dependence on the geometry 

In this section we outline the behaviour of the attractors as functions of one of the 
parameters which determine the shape of the billiard: s x . We take two electric fields (one 
large and one small) for which we had found apparently chaotic behaviours. 
Let e = 0.374 be the electric field, L = 1.291 and s y — 1.1; we study the attractor reliance on 
s x G [0, 1] with step 0.01. We considered trajectories of length 10 7 collisions and looked at the 
last 10 3 points. We found a complicated scenario, analogous to the one for the dependence 
on e, as s x was varied. In the chosen range, there are more cases with apparently chaotic 
attractors than with (quasi) periodic steady states. Increasing the number of collisions in 
each trajectory up to 5 x 10 7 , we find that apparent chaos persists, in other words it seems 
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that the effect of the electric field prevails on the geometry effects. To analyze more carefully 
this fact, we have honed the interval [0.09,0.1] for s x , with step 0.0005 and for trajectories 
of 10 8 collisions. A sudden transition between apparent chaos and periodic orbits happens 
at s x = 0.099. 

For e = 0.014, s x was taken in [0, 1] with step 0.01, for trajectories of length 10 7 collisions. 
We found that the situation is slightly different from the small field case. The apparently 
chaotic steady states occupy a smaller region of the phase space than the steady states of 
the e = 0.374 case. For instance, the asymptotic state of e = 0.014 seems to fill almost 
completely the space 6, differently from the case of e = 0.374. Increasing the length of 
the trajectories up to 5 x 10 7 collisions, various apparently chaotic cases reduce to (quasi) 
periodic orbits, while other survive. Finally, we increased the number of collisions up to 10 8 
for e = 0.014 again finding that some apparently chaotic cases turned into quasi periodic 
cases, while other survived. 

We conclude that the dependence of the attractors on the parameter s x is qualitatively 
similar to the dependence on e, although at times the dependence on e prevails. 



VI. CONCLUSION 



In this paper we have examined the non equilibrium version of the Ehrenfest gas, which 
is a billiard model with an electric field and a Gaussian thermostat, whose point particle 
moves in the plane and undergoes elastic collisions with rhomboidal obstacles. The moti- 
vation was to understand the dynamics of this nonequilibrium model, whose obstacles have 
flat surfaces, hence do not defocus the trajectories, while its thermostat, which makes dis- 
sipative the dynamics, does focus them j26f |. We have shown that periodic orbits with one 
positive Lyapunov exponent embedded in what appear to be chaotic attractors, do exist. 
Our numerical results have identified electric fields whose dynamics is strongly suggested 
to be chaotic, although conclusive results are out of reach at present, because of a very 
peculiar phenomena. The stationary state, even if attracting and trivial, requires very long 
times to be reached, because it coexists with a transiently chaotic state which covers most 
of the phase space. The dependence on the model parameters of the steady state behaviour 
also presents peculiar features which, to the best of our knowledge, have not been observed 
before: a very irregular, possibly discontinuous dependence of the attracting orbit, and/or 
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of the size and shape of the stable islands, on e and s x . The peculiarity remains even if it 
is eventually proved that all steady states are periodic, because of their coexistence with 
transiently chaotic states, and because of the consequent irregular behaviour of the conver- 
gence rates. As the sum of the Lyapunov exponents is proportional to the travelled distance, 
cf. Eq. (TTTj) . very irregular transport properties are obtained as well, perhaps as irregular as 
those conjectured for other low- dimensional dynamical systems 22]. 
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Figures 




FIG. 1: The modified Ehrenfest gas. In the present paper, the side of the elementary cell is set to 
1.291, while the semiaxis of the rhombus are chosen to be 1.1 and 0.7573 respectively. 
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FIG. 2: The closed periodic-two orbit SF 2 = (K-i _i, r 4 )(^i,i, r 3 ) for = ^, x; = (left panel). 
The open periodic orbit = (7^0,0, t"4)(^i,i, ^3)(^2,o, ^4) for = 2, = 0.8, = 0.5 (right 
panel). 
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FIG. 4: The last 10 4 iterates of the bounce map for e = 0.374 out of a trajectory of 5 x 10 7 collisions 
(left panel) and for a trajectory of 1.5 x 10 8 collisions, starting from the last phase space point 
of the previous trajectory (right panel). The distribution appears to be the same, suggesting that 
the system is in an apparently stationary state. Here, r\ = cos (p, with <p the angle between the 
outgoing velocity and the side of the rhombus, and r is the perimetral distance of the collision 
point from the right corner of the rhombus (77 and r are called Birkhoff coordinates). 
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FIG. 5: The open period-four orbit (7^o,o> r 3)C^-i,-i> r i)C^i,-i; r 4)(^-o,0) r 2)(7^2,0) r 3) m the plane 
for e = 0.374, s y = 0.7573, s x = 1.1, L = 1.291, 9 = -2.0619, x = -0.4885 (left panel). The four 
points of this periodic orbit, are inflated to show their embedding in the attractor (right panel). 
The coordinates are as in Figure 4. 
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FIG. 6: The last 7000 points of the bounce map for e = 0.5 out of a trajectory of length 5 x 10 
collisions (left panel) and the last 10 4 points of a trajectory of length 2 x 10 8 with same e (right 
panel). The coordinates are defined as in Figure 4. 



20 



0.01 - i \ ! 

All 

- * 







/" - 








^ II - 






















/ 













\ i 






— i — 

e 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 



FIG. 7: The largest Lyapunov exponent, Ai, for different lengths of trajectory: 10 7 (left panel) 
and 3 x 10 8 (right panel) for e G [0.002,0.1] and random initial condition. 
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FIG. 8: Behaviour of the finite time Lyapunov exponents Ai t 2 with the number of collisions n, for 
e = 0.087 (left panel). The logarithmic scale in n clearly separates the initial, intermediate and 
asymptotic regimes. The exponents rapidly converge to one positive and one negative value, which 
persist for very long, up to a number iV c of collisions, dependent on the initial condition. After 
N c , both exponents converge as 1/n to the value —0.004838. For random initial conditions, N c 
is of order O(10 8 ), for initial conditions close to the asymptotic periodic orbit (right panel), N c 
varies irregularly between O(10) and O(10 7 ). The right panel reports the last 10 4 points of 3 • 10 s 
collisions, which testifies that the motion has settled on a periodic orbit of 19 points only. The 
coordinates are defined as in Figure 4. 
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FIG. 9: The upper half of the projection of the multi-furcation diagram onto the 6 axis, for e in 
the range [0.01; 1]. The last 10 3 points of trajectories with length 10 7 collisions have been plotted 
(left panel). The multi-furcation diagram in the (e; 9) plane, for e in the range [1,1.3], from the 
last 10 3 collisions out of a trajectory of 10 7 collisions (right panel). 



Tables 



Field 


Collisions 


Lyapunov exponent 


0.374 


10 6 


0.144232 


0.374 


1.5 x 10 7 


0.144317 


0.374 


6.5 x 10 7 


0.144291 


0.374 


2.15 x 10 8 


0.144320 


0.5 


5 x 10 7 


0.166648 


0.5 


2 x 10 8 


0.166622 



TABLE I: Numerically computed Lyapunov exponents for different trajectories and fields. 



Field 


x 


Vo 


#0 


Ai 


A 2 


0.374 


-0.48971578282741385 


-0.3886737605834475 


-2.06199461833 


0.108316 


-0.281452 


0.5 


0.30576674801010 


0.6558650167554337 


-0.31873995693500 


0.155569 


-0.384674 



TABLE II: Initial conditions and Lyapunov exponents of periodic orbits of period four. All data 
have been computed analitically. 



